Cervezas Artesanas

1 Introducción

Debido a mi afición por la Cerveza Artesana, he seleccionado de Kaggle el dataset “beer-recipes” https://www.kaggle.com/jtrofe/beer-recipes

1.1 Descripción:

Este es un conjunto de datos de 75.000 cervezas caseras con más de 176 estilos diferentes. Los registros de cerveza son reportados por los usuarios y se clasifican de acuerdo a uno de los 176 estilos diferentes. Estas recetas entran tanto o tan poco en detalle como el usuario proporcionó, pero hay al menos 5 columnas útiles donde se introdujeron datos para cada una: Gravedad Original, Gravedad Final, ABV, IBU, y Color

1.2 Contenido

Datos básicos de las recetas de cerveza presentadas por los usuarios y raspadas de Brewer’s Friend. Todas las columnas están estandarizadas excepto “Método de cebado” y “Cantidad de cebado”, que parece que permite a los usuarios escribir lo que quieran. Agradecimientos

El sitio Brewer’s Friend permite a los usuarios compartir sus recetas de cerveza casera. Este conjunto de datos contiene una selección de las recetas subidas hasta ahora. Inspiración

1.3 ¿Qué hay en la cerveza casera?

Sería interesante ver si los datos proporcionados son suficientes para definir cada clase o si hay patrones no descubiertos. En el futuro podría ser posible revisar y raspar información más detallada para cada receta, como la levadura y los lúpulos específicos utilizados.

1.4 Detalle de las variables Principales:

  • ABV (Alcohol By Volume)

Alcohol por Volumen. Expresado como porcentaje (%), se utiliza como la medida de la intensidad del alcohol, basado en la proporción de su contenido por volumen de cerveza.

  • stilos de Cerveza

Los estilos de cerveza son las categorías mediante las cuales se identifican y clasifican las cervezas en base a sus características de apariencia, aroma, sabor y sensación en boca, establecidas mediante un rango aproximado de estadísticas vitales.

  • IBU (International Bitterness Unit)

Unidad Internacional de amargor. Es la medida utilizada usada para definir el amargor de una cerveza. Un IBU equivale a un miligramo de iso-alfa-ácidos disueltos en un litro de cerveza. Cuanto mayor es el valor, más amarga es la cerveza.

  • Color:

El color de la cerveza, depende principalmente, del tipo o tipos de maltas que se utilizan durante su elaboración. Hay 2 metodos el Americano (SRM) y el europeo (EBC)

2 El Objetivo:

Analizar si hay relación de los ingredientes y caracteristicas de cada cerveza que podamos obtener un patron de clasificación de las mismas, a pesar de su estilo.

3 Limpieza y preparación entorno trabajo

# Limpiamos el entorno de Trabajo
rm(list=ls())

# Limpiamos la consola
cat("\014")
# Comprobamos que está bien establecido el directorio
getwd()
## [1] "/home/oscar/Documentos/R/Cerveza"
dir()
## [1] "cerveza-cluster.Rmd" "data"                "header.html"
#indicamos el directorio de trabajo
setwd("~/Documentos/R/Cerveza")

4 Carga de librerias

packages <- c("psych", "dplyr", "skimr", "magrittr", "stringr", "epiDisplay", "summarytools","GGally", "ggplot2", "corrplot" , "githubinstall", "e1071", "purrr")
newpack  = packages[!(packages %in% installed.packages()[,"Package"])]

if(length(newpack)) install.packages(newpack)
a=lapply(packages, library, character.only=TRUE)

5 Descarga Database

Utilizaré un dataset proveniente de “kaggle”

Debido a la restricción que ofrece para poder descargar vía R su Dataset, voy a utilizar la terminal y mediante bash lo voy a descargar y descomprimir.

Este semestre estoy tambien matriculado en “Programación en Scripting”

Para obtener la url y descargar el dataset utilizo la extensión para Chrome wget/curl y realizo la siguientes etapas:

1 - Instalo la extensión.

2 - Empiezo a descargar y lo paro.

3 - Consigo el link de la extensión.

4 - wget o curl + el link.

Descargo el fichero: Instalar kaggle - utilizado el terminal mediante comandos bash.



# pip install kaggle

Obtengo el dataset - utilizado el terminal mediante comandos bash.

Eliminio los posibles ficheros que pudiera tener y descomprimo. - utilizado el terminal mediante comandos bash.

Descomprimo.

unzip(zipfile = "/home/oscar/Descargas/beer-recipes.zip", exdir = "/home/oscar/Documentos/R/Cerveza/data")

El dataset que voy a utilizar es el “recipeData.csv” de 13mb. No utilizo el styleData ya que dichos datos están incorporados en el anterior.

recipeData <- read.csv("~/Documentos/R/Cerveza/data/recipeData.csv", encoding = "latin1")
head(recipeData)
##   BeerID                            Name
## 1      1               Vanilla Cream Ale
## 2      2     Southern Tier Pumking clone
## 3      3     Zombie Dust Clone - EXTRACT
## 4      4   Zombie Dust Clone - ALL GRAIN
## 5      5 Bakke Brygg Belgisk Blonde 50 l
## 6      6    Sierra Nevada Pale Ale Clone
##                                                           URL
## 1                /homebrew/recipe/view/1633/vanilla-cream-ale
## 2     /homebrew/recipe/view/16367/southern-tier-pumking-clone
## 3        /homebrew/recipe/view/5920/zombie-dust-clone-extract
## 4      /homebrew/recipe/view/5916/zombie-dust-clone-all-grain
## 5 /homebrew/recipe/view/89534/bakke-brygg-belgisk-blonde-50-l
## 6    /homebrew/recipe/view/28546/sierra-nevada-pale-ale-clone
##                                Style StyleID Size.L.    OG    FG  ABV   IBU
## 1                          Cream Ale      45   21.77 1.055 1.013 5.48 17.65
## 2 Holiday/Winter Special Spiced Beer      85   20.82 1.083 1.021 8.16 60.65
## 3                       American IPA       7   18.93 1.063 1.018 5.91 59.25
## 4                       American IPA       7   22.71 1.061 1.017 5.80 54.48
## 5                  Belgian Blond Ale      20   50.00 1.060 1.010 6.48 17.84
## 6                  American Pale Ale      10   24.61 1.055 1.013 5.58 40.12
##   Color BoilSize BoilTime BoilGravity Efficiency MashThickness       SugarScale
## 1  4.83    28.39       75       1.038         70           N/A Specific Gravity
## 2 15.64    24.61       60        1.07         70           N/A Specific Gravity
## 3  8.98    22.71       60         N/A         70           N/A Specific Gravity
## 4  8.50    26.50       60         N/A         70           N/A Specific Gravity
## 5  4.57    60.00       90        1.05         72           N/A Specific Gravity
## 6  8.00    29.34       70       1.047         79           N/A Specific Gravity
##   BrewMethod PitchRate PrimaryTemp PrimingMethod  PrimingAmount UserId
## 1  All Grain       N/A       17.78    corn sugar         4.5 oz    116
## 2  All Grain       N/A         N/A           N/A            N/A    955
## 3    extract       N/A         N/A           N/A            N/A     NA
## 4  All Grain       N/A         N/A           N/A            N/A     NA
## 5  All Grain       N/A          19    Sukkerlake 6-7 g sukker/l  18325
## 6  All Grain         1         N/A           N/A            N/A   5889
print("Nombre de  las columnas: ")
## [1] "Nombre de  las columnas: "
names(recipeData)
##  [1] "BeerID"        "Name"          "URL"           "Style"        
##  [5] "StyleID"       "Size.L."       "OG"            "FG"           
##  [9] "ABV"           "IBU"           "Color"         "BoilSize"     
## [13] "BoilTime"      "BoilGravity"   "Efficiency"    "MashThickness"
## [17] "SugarScale"    "BrewMethod"    "PitchRate"     "PrimaryTemp"  
## [21] "PrimingMethod" "PrimingAmount" "UserId"
print("Clase: ")
## [1] "Clase: "
class(recipeData)
## [1] "data.frame"

6 Exploración

print("Dimensión: ")
## [1] "Dimensión: "
dim(recipeData)
## [1] 73814    23
print("Número de columnas: ")
## [1] "Número de columnas: "
length(recipeData)
## [1] 23
print("Estructura: ")
## [1] "Estructura: "
str(recipeData)
## 'data.frame':    73814 obs. of  23 variables:
##  $ BeerID       : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Name         : chr  "Vanilla Cream Ale" "Southern Tier Pumking clone" "Zombie Dust Clone - EXTRACT" "Zombie Dust Clone - ALL GRAIN" ...
##  $ URL          : chr  "/homebrew/recipe/view/1633/vanilla-cream-ale" "/homebrew/recipe/view/16367/southern-tier-pumking-clone" "/homebrew/recipe/view/5920/zombie-dust-clone-extract" "/homebrew/recipe/view/5916/zombie-dust-clone-all-grain" ...
##  $ Style        : chr  "Cream Ale" "Holiday/Winter Special Spiced Beer" "American IPA" "American IPA" ...
##  $ StyleID      : int  45 85 7 7 20 10 86 45 129 86 ...
##  $ Size.L.      : num  21.8 20.8 18.9 22.7 50 ...
##  $ OG           : num  1.05 1.08 1.06 1.06 1.06 ...
##  $ FG           : num  1.01 1.02 1.02 1.02 1.01 ...
##  $ ABV          : num  5.48 8.16 5.91 5.8 6.48 5.58 7.09 5.36 5.77 8.22 ...
##  $ IBU          : num  17.6 60.6 59.2 54.5 17.8 ...
##  $ Color        : num  4.83 15.64 8.98 8.5 4.57 ...
##  $ BoilSize     : num  28.4 24.6 22.7 26.5 60 ...
##  $ BoilTime     : int  75 60 60 60 90 70 90 75 75 60 ...
##  $ BoilGravity  : chr  "1.038" "1.07" "N/A" "N/A" ...
##  $ Efficiency   : num  70 70 70 70 72 79 75 70 73 70 ...
##  $ MashThickness: chr  "N/A" "N/A" "N/A" "N/A" ...
##  $ SugarScale   : chr  "Specific Gravity" "Specific Gravity" "Specific Gravity" "Specific Gravity" ...
##  $ BrewMethod   : chr  "All Grain" "All Grain" "extract" "All Grain" ...
##  $ PitchRate    : chr  "N/A" "N/A" "N/A" "N/A" ...
##  $ PrimaryTemp  : chr  "17.78" "N/A" "N/A" "N/A" ...
##  $ PrimingMethod: chr  "corn sugar" "N/A" "N/A" "N/A" ...
##  $ PrimingAmount: chr  "4.5 oz" "N/A" "N/A" "N/A" ...
##  $ UserId       : int  116 955 NA NA 18325 5889 1051 116 116 NA ...

Primera eliminación de columnas no necesarias

data <- dplyr::select(recipeData,-BeerID, -Name, -URL, -PrimingMethod, -PrimingMethod,-UserId, -StyleID)

summary(data)
##     Style              Size.L.              OG               FG        
##  Length:73814       Min.   :   1.00   Min.   : 1.000   Min.   :-0.003  
##  Class :character   1st Qu.:  18.93   1st Qu.: 1.051   1st Qu.: 1.011  
##  Mode  :character   Median :  20.82   Median : 1.058   Median : 1.013  
##                     Mean   :  43.94   Mean   : 1.406   Mean   : 1.076  
##                     3rd Qu.:  23.66   3rd Qu.: 1.069   3rd Qu.: 1.017  
##                     Max.   :9200.00   Max.   :34.035   Max.   :23.425  
##                     NA's   :1         NA's   :1        NA's   :1       
##       ABV              IBU              Color           BoilSize      
##  Min.   : 0.000   Min.   :   0.00   Min.   :  0.00   Min.   :   1.00  
##  1st Qu.: 5.080   1st Qu.:  23.37   1st Qu.:  5.17   1st Qu.:  20.82  
##  Median : 5.790   Median :  35.77   Median :  8.44   Median :  27.44  
##  Mean   : 6.137   Mean   :  44.28   Mean   : 13.41   Mean   :  49.73  
##  3rd Qu.: 6.830   3rd Qu.:  56.39   3rd Qu.: 16.79   3rd Qu.:  30.00  
##  Max.   :54.720   Max.   :3409.30   Max.   :186.00   Max.   :9700.00  
##  NA's   :1        NA's   :1         NA's   :1        NA's   :1        
##     BoilTime      BoilGravity          Efficiency     MashThickness     
##  Min.   :  0.00   Length:73814       Min.   :  0.00   Length:73814      
##  1st Qu.: 60.00   Class :character   1st Qu.: 65.00   Class :character  
##  Median : 60.00   Mode  :character   Median : 70.00   Mode  :character  
##  Mean   : 65.08                      Mean   : 66.35                     
##  3rd Qu.: 60.00                      3rd Qu.: 75.00                     
##  Max.   :240.00                      Max.   :100.00                     
##  NA's   :1                           NA's   :1                          
##   SugarScale         BrewMethod         PitchRate         PrimaryTemp       
##  Length:73814       Length:73814       Length:73814       Length:73814      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  PrimingAmount     
##  Length:73814      
##  Class :character  
##  Mode  :character  
##                    
##                    
##                    
## 
glimpse(data)
## Rows: 73,814
## Columns: 17
## $ Style         <chr> "Cream Ale", "Holiday/Winter Special Spiced Beer", "Ame…
## $ Size.L.       <dbl> 21.77, 20.82, 18.93, 22.71, 50.00, 24.61, 22.71, 20.82,…
## $ OG            <dbl> 1.055, 1.083, 1.063, 1.061, 1.060, 1.055, 1.072, 1.054,…
## $ FG            <dbl> 1.013, 1.021, 1.018, 1.017, 1.010, 1.013, 1.018, 1.014,…
## $ ABV           <dbl> 5.48, 8.16, 5.91, 5.80, 6.48, 5.58, 7.09, 5.36, 5.77, 8…
## $ IBU           <dbl> 17.65, 60.65, 59.25, 54.48, 17.84, 40.12, 268.71, 19.97…
## $ Color         <dbl> 4.83, 15.64, 8.98, 8.50, 4.57, 8.00, 6.33, 5.94, 34.76,…
## $ BoilSize      <dbl> 28.39, 24.61, 22.71, 26.50, 60.00, 29.34, 30.28, 28.39,…
## $ BoilTime      <int> 75, 60, 60, 60, 90, 70, 90, 75, 75, 60, 90, 90, 60, 60,…
## $ BoilGravity   <chr> "1.038", "1.07", "N/A", "N/A", "1.05", "1.047", "N/A", …
## $ Efficiency    <dbl> 70, 70, 70, 70, 72, 79, 75, 70, 73, 70, 74, 70, 70, 30,…
## $ MashThickness <chr> "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "1.4",…
## $ SugarScale    <chr> "Specific Gravity", "Specific Gravity", "Specific Gravi…
## $ BrewMethod    <chr> "All Grain", "All Grain", "extract", "All Grain", "All …
## $ PitchRate     <chr> "N/A", "N/A", "N/A", "N/A", "N/A", "1", "N/A", "N/A", "…
## $ PrimaryTemp   <chr> "17.78", "N/A", "N/A", "N/A", "19", "N/A", "N/A", "N/A"…
## $ PrimingAmount <chr> "4.5 oz", "N/A", "N/A", "N/A", "6-7 g sukker/l", "N/A",…
# library(psych)
describe(data)
##                vars     n  mean     sd median trimmed   mad    min     max
## Style*            1 73814   NaN     NA     NA     NaN    NA    Inf    -Inf
## Size.L.           2 73813 43.94 180.43  20.82   22.33  2.80   1.00 9200.00
## OG                3 73813  1.41   2.20   1.06    1.06  0.01   1.00   34.03
## FG                4 73813  1.08   0.43   1.01    1.01  0.00   0.00   23.42
## ABV               5 73813  6.14   1.88   5.79    5.94  1.22   0.00   54.72
## IBU               6 73813 44.28  42.95  35.77   39.31 21.93   0.00 3409.30
## Color             7 73813 13.41  11.95   8.44   11.11  6.17   0.00  186.00
## BoilSize          8 73813 49.73 193.31  27.44   26.63  6.76   1.00 9700.00
## BoilTime          9 73813 65.08  15.02  60.00   63.37  0.00   0.00  240.00
## BoilGravity*     10 73814  1.35   1.93   1.05    1.05  0.01   0.00   52.60
## Efficiency       11 73813 66.35  14.09  70.00   68.43  7.41   0.00  100.00
## MashThickness*   12 73814  2.13   1.68   1.50    1.95  0.37   0.00  100.00
## SugarScale*      13 73814   NaN     NA     NA     NaN    NA    Inf    -Inf
## BrewMethod*      14 73814   NaN     NA     NA     NaN    NA    Inf    -Inf
## PitchRate*       15 73814  0.75   0.39   0.75    0.70  0.37   0.00    2.00
## PrimaryTemp*     16 73814 19.18   4.22  20.00   19.39  1.65 -17.78  114.00
## PrimingAmount*   17 73813 64.33 133.36   8.00   46.65 10.75   0.00 2500.00
##                  range  skew kurtosis   se
## Style*            -Inf    NA       NA   NA
## Size.L.        9199.00 15.87   397.42 0.66
## OG               33.03  6.70    47.28 0.01
## FG               23.43  9.78   174.25 0.00
## ABV              54.72  4.95    86.81 0.01
## IBU            3409.30 19.20  1113.91 0.16
## Color           186.00  1.59     2.20 0.04
## BoilSize       9699.00 15.85   402.07 0.71
## BoilTime        240.00  1.38    12.60 0.06
## BoilGravity*     52.60  7.38    67.09 0.01
## Efficiency      100.00 -1.48     1.79 0.05
## MashThickness*  100.00 16.85   540.01 0.01
## SugarScale*       -Inf    NA       NA   NA
## BrewMethod*       -Inf    NA       NA   NA
## PitchRate*        2.00  1.05     0.74 0.00
## PrimaryTemp*    131.78  3.26    63.79 0.02
## PrimingAmount* 2500.00 11.41   190.89 0.49
# valores vacios
print("Mostrar variables con campos na")
## [1] "Mostrar variables con campos na"
which (is.na(data))
##  [1]  147628  221442  295256  369070  442884  516698  590512  664326  811954
## [10] 1182422
print("Es cierto que hay valores  na?")
## [1] "Es cierto que hay valores  na?"
any(is.na(data))
## [1] TRUE
print("Suma valores na")
## [1] "Suma valores na"
sum(is.na(data))
## [1] 10
print("Mostrar variables con campos na")
## [1] "Mostrar variables con campos na"
colSums(is.na(data))
##         Style       Size.L.            OG            FG           ABV 
##             0             1             1             1             1 
##           IBU         Color      BoilSize      BoilTime   BoilGravity 
##             1             1             1             1             0 
##    Efficiency MashThickness    SugarScale    BrewMethod     PitchRate 
##             1             0             0             0             0 
##   PrimaryTemp PrimingAmount 
##             0             1
print("Mostrar variables con datos N/A")
## [1] "Mostrar variables con datos N/A"
colSums(data=="N/A")
##         Style       Size.L.            OG            FG           ABV 
##           595            NA            NA            NA            NA 
##           IBU         Color      BoilSize      BoilTime   BoilGravity 
##            NA            NA            NA            NA          2990 
##    Efficiency MashThickness    SugarScale    BrewMethod     PitchRate 
##            NA         29843             0             0         39228 
##   PrimaryTemp PrimingAmount 
##         22646            NA
print("Valores con integrogación")
## [1] "Valores con integrogación"
colSums(data=="?")
##         Style       Size.L.            OG            FG           ABV 
##             0            NA            NA            NA            NA 
##           IBU         Color      BoilSize      BoilTime   BoilGravity 
##            NA            NA            NA            NA             0 
##    Efficiency MashThickness    SugarScale    BrewMethod     PitchRate 
##            NA             0             0             0             0 
##   PrimaryTemp PrimingAmount 
##             0            NA
print("Valores con espacios en blanco")
## [1] "Valores con espacios en blanco"
colSums(data==" ")
##         Style       Size.L.            OG            FG           ABV 
##             0            NA            NA            NA            NA 
##           IBU         Color      BoilSize      BoilTime   BoilGravity 
##            NA            NA            NA            NA             0 
##    Efficiency MashThickness    SugarScale    BrewMethod     PitchRate 
##            NA             0             0             0             0 
##   PrimaryTemp PrimingAmount 
##             0            NA
data$Style <- gsub("N/A", NA, data$Style, fixed = TRUE)
data$BoilSize <- gsub("N/A", NA, data$BoilSize, fixed = TRUE)
data$PitchRate <- gsub("N/A",NA, data$PitchRate, fixed = TRUE)
data$Style <- gsub("N/A",NA, data$Style, fixed = TRUE)
data$BoilGravity <- gsub("N/A",NA, data$BoilGravity, fixed = TRUE)
# me quedo con los datos completos

completos <- complete.cases(df)

datos <- df[completos,]
# valores vacios
print("Mostrar variables con campos na")
## [1] "Mostrar variables con campos na"
which (is.na(datos))
## integer(0)
print("Es cierto que hay valores  na?")
## [1] "Es cierto que hay valores  na?"
any(is.na(datos))
## [1] FALSE
print("Suma valores na")
## [1] "Suma valores na"
sum(is.na(datos))
## [1] 0
print("Mostrar variables con campos na")
## [1] "Mostrar variables con campos na"
colSums(is.na(datos))
##       Style     Size.L.          OG          FG         ABV         IBU 
##           0           0           0           0           0           0 
##       Color    BoilSize    BoilTime BoilGravity  Efficiency  SugarScale 
##           0           0           0           0           0           0 
##  BrewMethod 
##           0
print("Mostrar variables con datos N/A")
## [1] "Mostrar variables con datos N/A"
colSums(datos=="N/A")
##       Style     Size.L.          OG          FG         ABV         IBU 
##           0           0           0           0           0           0 
##       Color    BoilSize    BoilTime BoilGravity  Efficiency  SugarScale 
##           0           0           0           0           0           0 
##  BrewMethod 
##           0
print("Valores con integrogación")
## [1] "Valores con integrogación"
colSums(datos=="?")
##       Style     Size.L.          OG          FG         ABV         IBU 
##           0           0           0           0           0           0 
##       Color    BoilSize    BoilTime BoilGravity  Efficiency  SugarScale 
##           0           0           0           0           0           0 
##  BrewMethod 
##           0
print("Valores con espacios en blanco")
## [1] "Valores con espacios en blanco"
colSums(datos==" ")
##       Style     Size.L.          OG          FG         ABV         IBU 
##           0           0           0           0           0           0 
##       Color    BoilSize    BoilTime BoilGravity  Efficiency  SugarScale 
##           0           0           0           0           0           0 
##  BrewMethod 
##           0
# Modificar los estilos en mayuscula
# library(stringr)
datos$Style <- str_to_upper(datos$Style, locale = "en")  #vocablos en ingles
# el resto de vairables character en minusculas
datos$BoilSize <- str_to_lower(datos$BoilSize, locale = "en")
datos$BoilGravity <- str_to_lower(datos$BoilGravity, locale = "en")
# Elimino espacios en blanco
datos$Style <- str_trim(datos$Style, side = "both")
# Separo la columna Style en 2 partes

datos <- datos %>%
  tidyr::separate(Style, c("Estilo", "Clase", "Otros"), sep = " ", fill="right")

# Elimino la columna "Otros"
datos <- dplyr::select(datos, -Otros)
# Sustituyo valores na de la columna clase con los valores de la columna Estilo
datos$Clase[is.na(datos$Clase)]<- as.character(datos$Estilo[is.na(datos$Clase)])
# Ref. https://stackoverflow.com/questions/15629885/replace-na-in-column-with-value-in-adjacent-column

# Relleno el espacio en blanco en Clase indicando "otros"
datos$Clase <- sub("^$", "OTHERS", datos$Clase)

# Ref. https://stackoverflow.com/questions/21243588/replace-blank-cells-with-character

# Elimino caracteres especiales

datos$Clase <- gsub("[[:punct:]]", "", datos$Clase)

Reducción dataset en 5 tipos de cerverza

# Seleccionamos las Clases de cerveza por nº de frecuencia, hasta que llega a 50% del total y poder obtener clusteres bien definicos
top_data1 <- filter(datos, Clase %in% c("IPA","PALE","STOUT","ALE","PORTER"))

Reduccion de observaciones

# Dado que el dataset es de tamaño muy grande, reduzco a 1000 primeras filas aleatorias sin reposición


top_data <- top_data1[sample(nrow(top_data1), 1000, replace=FALSE),]

Nuevo resumen de datos

skimr::skim(top_data)
Data summary
Name top_data
Number of rows 1000
Number of columns 14
_______________________
Column type frequency:
character 6
numeric 8
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
Estilo 0 1 3 13 0 20 0
Clase 0 1 3 6 0 5 0
BoilSize 0 1 1 7 0 207 0
BoilGravity 0 1 2 5 0 135 0
SugarScale 0 1 5 16 0 2 0
BrewMethod 0 1 4 12 0 4 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Size.L. 0 1 36.81 116.43 3.03 18.93 20.82 23.00 1994.91 ▇▁▁▁▁
OG 0 1 1.30 1.76 1.02 1.05 1.06 1.07 16.54 ▇▁▁▁▁
FG 0 1 1.06 0.36 1.00 1.01 1.01 1.02 6.48 ▇▁▁▁▁
ABV 0 1 6.23 1.42 2.25 5.26 6.05 6.93 15.92 ▂▇▁▁▁
IBU 0 1 60.29 62.38 0.00 33.94 49.63 72.28 1605.83 ▇▁▁▁▁
Color 0 1 14.07 12.98 2.25 5.87 7.92 14.65 50.00 ▇▁▁▁▁
BoilTime 0 1 64.51 13.09 20.00 60.00 60.00 60.00 240.00 ▇▂▁▁▁
Efficiency 0 1 66.45 13.80 25.00 65.00 70.00 75.00 95.00 ▂▁▂▇▁

Convierto atrubutos character a factor y a numeric

# convierto la columna character a factor
top_data$Estilo <- as.factor(top_data$Estilo)
top_data$Clase <- as.factor(top_data$Clase)
top_data$BoilSize <- as.factor(top_data$BoilSize)
top_data$BoilGravity <- as.factor(top_data$BoilGravity)
# convierto variables con mucho números únicos en númericos
top_data$BoilSize <- as.numeric(top_data$BoilSize)
top_data$BoilGravity <- as.numeric(top_data$BoilGravity)

Busco valores únicos

unique(top_data$Estilo)
##  [1] AMERICAN      DRY           BALTIC        IMPERIAL      SWEET        
##  [6] BELGIAN       BLONDE        CREAM         SPECIALTY     ROBUST       
## [11] CZECH         ENGLISH       BROWN         KELLERBIER:   OATMEAL      
## [16] DOUBLE        OLD           IRISH         INTERNATIONAL TROPICAL     
## 20 Levels: AMERICAN BALTIC BELGIAN BLONDE BROWN CREAM CZECH DOUBLE ... TROPICAL
unique(top_data$Clase)
## [1] IPA    PALE   STOUT  PORTER ALE   
## Levels: ALE IPA PALE PORTER STOUT
length(unique(unlist(top_data[c("Estilo")])))
## [1] 20
length(unique(unlist(top_data[c("Clase")])))
## [1] 5
agg <- aggregate(Estilo ~ Clase, top_data, function(x) length(unique(x)))
agg
##    Clase Estilo
## 1    ALE      3
## 2    IPA      5
## 3   PALE      5
## 4 PORTER      5
## 5  STOUT      7

Matriz

pairs.panels(top_data, pch=21,main="Gráfico 01.6: Matriz de Dispersión, Histograma y Correlación")

summarytools::freq(top_data$Clase, order = "freq")
## Frequencies  
## top_data$Clase  
## Type: Factor  
## 
##                Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ------------ ------ --------- -------------- --------- --------------
##          IPA    495     49.50          49.50     49.50          49.50
##         PALE    241     24.10          73.60     24.10          73.60
##        STOUT    122     12.20          85.80     12.20          85.80
##          ALE     76      7.60          93.40      7.60          93.40
##       PORTER     66      6.60         100.00      6.60         100.00
##         <NA>      0                               0.00         100.00
##        Total   1000    100.00         100.00    100.00         100.00
# library(epiDisplay)

epiDisplay::summ(top_data)
## 
## No. of observations = 1000
## 
##    Var. name   obs. mean   median  s.d.   min.   max.   
## 1  Estilo      1000 5.189  1       6.018  1      20     
## 2  Clase       1000 2.663  2       1.114  1      5      
## 3  Size.L.     1000 36.81  20.82   116.43 3.03   1994.91
## 4  OG          1000 1.3    1.06    1.76   1.02   16.54  
## 5  FG          1000 1.06   1.01    0.36   1      6.48   
## 6  ABV         1000 6.23   6.05    1.42   2.25   15.92  
## 7  IBU         1000 60.29  49.63   62.38  0      1605.83
## 8  Color       1000 14.07  7.92    12.98  2.25   50     
## 9  BoilSize    1000 91.19  95      47.34  1      207    
## 10 BoilTime    1000 64.51  60      13.09  20     240    
## 11 BoilGravity 1000 38.88  32.5    22.97  1      135    
## 12 Efficiency  1000 66.45  70      13.8   25     95     
## 13 SugarScale                                           
## 14 BrewMethod
tab1(top_data$Clase, sort.group = "decreasing", cum.percent = TRUE)

## top_data$Clase : 
##         Frequency Percent Cum. percent
## IPA           495    49.5         49.5
## PALE          241    24.1         73.6
## STOUT         122    12.2         85.8
## ALE            76     7.6         93.4
## PORTER         66     6.6        100.0
##   Total      1000   100.0        100.0
tab1(top_data$Estilo, sort.group = "decreasing", cum.percent = TRUE)

## top_data$Estilo : 
##               Frequency Percent Cum. percent
## AMERICAN            579    57.9         57.9
## SPECIALTY            65     6.5         64.4
## IMPERIAL             65     6.5         70.9
## BLONDE               48     4.8         75.7
## ENGLISH              34     3.4         79.1
## OATMEAL              27     2.7         81.8
## DOUBLE               27     2.7         84.5
## SWEET                25     2.5         87.0
## BELGIAN              23     2.3         89.3
## CREAM                21     2.1         91.4
## DRY                  18     1.8         93.2
## ROBUST               15     1.5         94.7
## INTERNATIONAL        13     1.3         96.0
## BROWN                10     1.0         97.0
## BALTIC                9     0.9         97.9
## CZECH                 8     0.8         98.7
## OLD                   7     0.7         99.4
## IRISH                 3     0.3         99.7
## TROPICAL              2     0.2         99.9
## KELLERBIER:           1     0.1        100.0
##   Total            1000   100.0        100.0
tab1(top_data$SugarScale, sort.group = "drecreasing", cum.percent = TRUE)

## top_data$SugarScale : 
##                  Frequency Percent Cum. percent
## Plato                   18     1.8          1.8
## Specific Gravity       982    98.2        100.0
##   Total               1000   100.0        100.0
tab1(top_data$BrewMethod, sort.group = "drecreasing", cum.percent = TRUE)

## top_data$BrewMethod : 
##              Frequency Percent Cum. percent
## All Grain          656    65.6         65.6
## BIAB               183    18.3         83.9
## extract            120    12.0         95.9
## Partial Mash        41     4.1        100.0
##   Total           1000   100.0        100.0

Guardo el Datamatrix

write.csv(top_data, file  =  "data/top_data.csv", row.names = FALSE)

Convierto el dataframe en Datamatrix

# Utilizando la libreria skimr obtenemos información adicional
dfmatrix <- top_data

cols <- c("Estilo","Clase","BoilSize","BoilGravity","SugarScale","BrewMethod")
# Convertimos lo valores character a factor
for (i in cols){
  dfmatrix[,i] <- as.numeric(dfmatrix[,i])
}

skim(dfmatrix)
Data summary
Name dfmatrix
Number of rows 1000
Number of columns 14
_______________________
Column type frequency:
numeric 14
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Estilo 0 1 5.19 6.02 1.00 1.00 1.00 10.00 20.00 ▇▁▂▁▁
Clase 0 1 2.66 1.11 1.00 2.00 2.00 3.00 5.00 ▁▇▃▁▂
Size.L. 0 1 36.81 116.43 3.03 18.93 20.82 23.00 1994.91 ▇▁▁▁▁
OG 0 1 1.30 1.76 1.02 1.05 1.06 1.07 16.54 ▇▁▁▁▁
FG 0 1 1.06 0.36 1.00 1.01 1.01 1.02 6.48 ▇▁▁▁▁
ABV 0 1 6.23 1.42 2.25 5.26 6.05 6.93 15.92 ▂▇▁▁▁
IBU 0 1 60.29 62.38 0.00 33.94 49.63 72.28 1605.83 ▇▁▁▁▁
Color 0 1 14.07 12.98 2.25 5.87 7.92 14.65 50.00 ▇▁▁▁▁
BoilSize 0 1 91.19 47.34 1.00 68.00 95.00 108.00 207.00 ▃▆▇▂▂
BoilTime 0 1 64.51 13.09 20.00 60.00 60.00 60.00 240.00 ▇▂▁▁▁
BoilGravity 0 1 38.88 22.97 1.00 26.00 32.50 43.00 135.00 ▅▇▁▁▁
Efficiency 0 1 66.45 13.80 25.00 65.00 70.00 75.00 95.00 ▂▁▂▇▁
SugarScale 1000 0 NaN NA NA NA NA NA NA
BrewMethod 1000 0 NaN NA NA NA NA NA NA
ggpairs(dfmatrix)

# Relación entre variables

correl <- cor(dfmatrix, method = "pearson")
print("Matriz de correlacion de Pearson")
## [1] "Matriz de correlacion de Pearson"
correl
##                   Estilo        Clase     Size.L.           OG          FG
## Estilo       1.000000000  0.237557775 -0.01088419  0.012448984  0.02799722
## Clase        0.237557775  1.000000000  0.01192920  0.026370989  0.05412004
## Size.L.     -0.010884189  0.011929199  1.00000000  0.127756217  0.09608860
## OG           0.012448984  0.026370989  0.12775622  1.000000000  0.93286139
## FG           0.027997220  0.054120040  0.09608860  0.932861393  1.00000000
## ABV          0.199401601 -0.003972842 -0.03665281 -0.024556935 -0.02974529
## IBU         -0.005751447 -0.136313374 -0.01806301 -0.023188529 -0.02644654
## Color        0.415466693  0.708621343 -0.01167523  0.009396223  0.03070150
## BoilSize     0.055167831 -0.029795743  0.05264564  0.026152204  0.03293212
## BoilTime     0.060180220  0.046733561  0.01739416  0.038584816  0.07311932
## BoilGravity  0.128738983  0.098392378  0.07629586  0.519083528  0.47878075
## Efficiency  -0.020434482 -0.026567755  0.15407467  0.103212289  0.09921376
## SugarScale            NA           NA          NA           NA          NA
## BrewMethod            NA           NA          NA           NA          NA
##                      ABV          IBU        Color     BoilSize    BoilTime
## Estilo       0.199401601 -0.005751447  0.415466693  0.055167831  0.06018022
## Clase       -0.003972842 -0.136313374  0.708621343 -0.029795743  0.04673356
## Size.L.     -0.036652813 -0.018063012 -0.011675228  0.052645638  0.01739416
## OG          -0.024556935 -0.023188529  0.009396223  0.026152204  0.03858482
## FG          -0.029745286 -0.026446541  0.030701499  0.032932118  0.07311932
## ABV          1.000000000  0.258307211  0.263626415  0.023477048  0.21078223
## IBU          0.258307211  1.000000000 -0.024002402 -0.014568568  0.07074333
## Color        0.263626415 -0.024002402  1.000000000  0.008327619  0.07208234
## BoilSize     0.023477048 -0.014568568  0.008327619  1.000000000  0.09745090
## BoilTime     0.210782230  0.070743333  0.072082343  0.097450898  1.00000000
## BoilGravity  0.369819088  0.083403195  0.194743035 -0.126576850 -0.00856133
## Efficiency   0.096397323 -0.018503861 -0.025605275  0.171964409  0.17102686
## SugarScale            NA           NA           NA           NA          NA
## BrewMethod            NA           NA           NA           NA          NA
##             BoilGravity  Efficiency SugarScale BrewMethod
## Estilo       0.12873898 -0.02043448         NA         NA
## Clase        0.09839238 -0.02656776         NA         NA
## Size.L.      0.07629586  0.15407467         NA         NA
## OG           0.51908353  0.10321229         NA         NA
## FG           0.47878075  0.09921376         NA         NA
## ABV          0.36981909  0.09639732         NA         NA
## IBU          0.08340320 -0.01850386         NA         NA
## Color        0.19474304 -0.02560527         NA         NA
## BoilSize    -0.12657685  0.17196441         NA         NA
## BoilTime    -0.00856133  0.17102686         NA         NA
## BoilGravity  1.00000000 -0.23641183         NA         NA
## Efficiency  -0.23641183  1.00000000         NA         NA
## SugarScale           NA          NA          1         NA
## BrewMethod           NA          NA         NA          1
# Buscamos  las variables con mayor correlación, hemos realizado una selección entre 0.40 y 0.85
correl_df <- as.data.frame(as.table(correl))
correl_df_a <- subset(correl_df, abs(Freq) > 0.50 & abs(Freq) <0.85)
  
correl_df_a[order(correl_df_a$Freq),]
##            Var1        Var2      Freq
## 53  BoilGravity          OG 0.5190835
## 144          OG BoilGravity 0.5190835
## 22        Color       Clase 0.7086213
## 100       Clase       Color 0.7086213
# Correlación cla
cor_clase <- cor(dfmatrix[-2], dfmatrix$Clase)
print(cor_clase)
##                     [,1]
## Estilo       0.237557775
## Size.L.      0.011929199
## OG           0.026370989
## FG           0.054120040
## ABV         -0.003972842
## IBU         -0.136313374
## Color        0.708621343
## BoilSize    -0.029795743
## BoilTime     0.046733561
## BoilGravity  0.098392378
## Efficiency  -0.026567755
## SugarScale            NA
## BrewMethod            NA
cor_clase1 <- as.data.frame(cor_clase)
cor_clase1[order(abs(cor_clase1$V1)),]
##  [1] -0.003972842  0.011929199  0.026370989 -0.026567755 -0.029795743
##  [6]  0.046733561  0.054120040  0.098392378 -0.136313374  0.237557775
## [11]  0.708621343           NA           NA

Salvo la relación con Estilo, seleccionolas que me parecen más importantes

Selecciono los mejores valores y reduzco el nº me observaciones

# Elegimos las mejores opciones e incorporo IBU y ABV

top_top <- dplyr::select(top_data, Clase, Color, BoilGravity, FG,  Efficiency, IBU, ABV)

Eliminamos los dataset que no necesitamos

rm(a, agg, data, datos, df, dfmatrix, recipeData, top_data1)

7 Outliers

Vamos detectar Outliers

par(mfrow=c(2,3))
boxplot(top_top$Color~top_top$Clase)
boxplot(top_top$BoilGravity~top_top$Clase)
boxplot(top_top$FG~top_top$Clase)
boxplot(top_top$Efficiency~top_top$Clase)
boxplot(top_top$IBU~top_top$Clase)
boxplot(top_top$ABV~top_top$Clase)

boxplot(ABV ~ Clase, 
        data = top_top,
        main = "ABV por Clase")

par(mfrow=c(3,2))
boxplot(top_top$Color~top_top$Clase)
boxplot(top_top$BoilGravity~top_top$Clase)
boxplot(top_top$FG~top_top$Clase)
boxplot(top_top$Efficiency~top_top$Clase)
boxplot(top_top$IBU~top_top$Clase)
boxplot(top_top$ABV~top_top$Clase)

# df es el dataFrame que recibimos (ej. activity)
# colNameData es la columna de los datos (ej. "steps")
# colNameBy es la columna por la que trocearemos (ej. "userId")
outliersReplace <- function(df, colNameData, colNameBy){
  # creamos una nueva columna llamada igual que colNameData pero con .R
  colNameData.R <- paste(colNameData, "R", sep=".")
  df[colNameData.R] <- df[colNameData]
  
  # obtenemos los IDs por los que partir el dataframe
  IDs <- unique(df[,c(colNameBy)])
  for (id in IDs){
    data <- df[df[colNameBy] == id, c(colNameData) ]
    
    Q  <- quantile(data)
    minimo <- Q[1]    # valor minimo
    Q1     <- Q[2]    # primer cuartil
    Me     <- Q[3]    # mediana
    Q3     <- Q[4]    # tercer cuartil
    maximo <- Q[5]    # valor maximo
    IQR    <- Q3 - Q1
    
    lowLimit  <- max(minimo, Q1 - 1.7*IQR)
    highLimit <- min(maximo, Q3 + 1.7*IQR)
    
    # todos los valores donde colNameBy es igual a id
    # y el valor de colNameData es > Q3 + 1.5 * IQR
    # lo reemplazamos por la mediana
    df[df[colNameBy] == id & df[colNameData] > highLimit, c(colNameData.R)] <- Me
    
    # lo mismo para el umbral inferior
    df[df[colNameBy] == id & df[colNameData] < lowLimit, c(colNameData.R)] <- Me
    
    cat(paste("El", colNameBy, id, "la mediana(", colNameData, ") ==", Me, "\n", sep=" " ))
    
  }
  df   # devolvemos el valor del dataFrame
}
top_top <- outliersReplace(top_top,"Color","Clase")
## El Clase IPA la mediana( Color ) == 7.66 
## El Clase PALE la mediana( Color ) == 6.47 
## El Clase STOUT la mediana( Color ) == 38.09 
## El Clase PORTER la mediana( Color ) == 33.185 
## El Clase ALE la mediana( Color ) == 4.635
top_top <- outliersReplace(top_top,"BoilGravity","Clase")
## El Clase IPA la mediana( BoilGravity ) == 35 
## El Clase PALE la mediana( BoilGravity ) == 28 
## El Clase STOUT la mediana( BoilGravity ) == 38 
## El Clase PORTER la mediana( BoilGravity ) == 30.5 
## El Clase ALE la mediana( BoilGravity ) == 25.5
top_top <- outliersReplace(top_top,"FG","Clase")
## El Clase IPA la mediana( FG ) == 1.014 
## El Clase PALE la mediana( FG ) == 1.011 
## El Clase STOUT la mediana( FG ) == 1.017 
## El Clase PORTER la mediana( FG ) == 1.0155 
## El Clase ALE la mediana( FG ) == 1.011
top_top <- outliersReplace(top_top,"Efficiency","Clase")
## El Clase IPA la mediana( Efficiency ) == 70 
## El Clase PALE la mediana( Efficiency ) == 70 
## El Clase STOUT la mediana( Efficiency ) == 70 
## El Clase PORTER la mediana( Efficiency ) == 70 
## El Clase ALE la mediana( Efficiency ) == 70
top_top <- outliersReplace(top_top,"IBU","Clase")
## El Clase IPA la mediana( IBU ) == 68.24 
## El Clase PALE la mediana( IBU ) == 39.39 
## El Clase STOUT la mediana( IBU ) == 37.9 
## El Clase PORTER la mediana( IBU ) == 33.62 
## El Clase ALE la mediana( IBU ) == 22.5
top_top <- outliersReplace(top_top,"ABV","Clase")
## El Clase IPA la mediana( ABV ) == 6.57 
## El Clase PALE la mediana( ABV ) == 5.34 
## El Clase STOUT la mediana( ABV ) == 5.845 
## El Clase PORTER la mediana( ABV ) == 6.02 
## El Clase ALE la mediana( ABV ) == 5.135

Vemos el resultado

par(mfrow = c(2,1))    # para ponerlos uno encima de otro
 
boxplot(ABV   ~ Clase, data = top_top, main = "Sin reemplazo")
boxplot(ABV.R ~ Clase, data = top_top, main = "Con reemplazo")

par(mfrow=c(3,2))
boxplot(top_top$Color.R~top_top$Clase)
boxplot(top_top$BoilGravity.R~top_top$Clase)
boxplot(top_top$FG.R~top_top$Clase)
boxplot(top_top$Efficiency.R~top_top$Clase)
boxplot(top_top$IBU.R~top_top$Clase)
boxplot(top_top$ABV.R~top_top$Clase)

head(top_top)
##       Clase Color BoilGravity      FG Efficiency   IBU  ABV Color.R
## 7608    IPA 16.72          26 1.01000         70 40.91 5.34   16.72
## 14197  PALE  5.09          32 1.00900         78 44.45 5.12    5.09
## 2761  STOUT 36.42         129 5.33594         66 44.52 4.88   36.42
## 29583   IPA  5.60          37 1.01500         70 71.78 6.27    5.60
## 1880  STOUT 50.00          25 1.00900         50 41.63 4.88   50.00
## 17415  PALE  9.43          31 1.01100         78 38.89 5.82    9.43
##       BoilGravity.R  FG.R Efficiency.R IBU.R ABV.R
## 7608             26 1.010           70 40.91  5.34
## 14197            32 1.009           78 44.45  5.12
## 2761             38 1.017           66 44.52  4.88
## 29583            37 1.015           70 71.78  6.27
## 1880             25 1.009           50 41.63  4.88
## 17415            31 1.011           78 38.89  5.82
top_top <- dplyr::select(top_top, -Color, -BoilGravity, -FG, -Efficiency, -IBU, -ABV)
names(top_top) <- c("Clase", "Color", "BoilGravity", "FG", "Efficiency", "IBU", "ABV")
str(top_top)
## 'data.frame':    1000 obs. of  7 variables:
##  $ Clase      : Factor w/ 5 levels "ALE","IPA","PALE",..: 2 3 5 2 5 3 4 2 2 4 ...
##  $ Color      : num  16.72 5.09 36.42 5.6 50 ...
##  $ BoilGravity: num  26 32 38 37 25 31 28 36 34 40 ...
##  $ FG         : num  1.01 1.01 1.02 1.01 1.01 ...
##  $ Efficiency : num  70 78 66 70 50 78 75 68 80 50 ...
##  $ IBU        : num  40.9 44.5 44.5 71.8 41.6 ...
##  $ ABV        : num  5.34 5.12 4.88 6.27 4.88 5.82 6.98 6.2 6.12 6.56 ...

8 Estructura de los datos transformados

skimr::skim(top_top)
Data summary
Name top_top
Number of rows 1000
Number of columns 7
_______________________
Column type frequency:
factor 1
numeric 6
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Clase 0 1 FALSE 5 IPA: 495, PAL: 241, STO: 122, ALE: 76

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Color 0 1 12.91 12.29 2.25 5.83 7.66 11.75 50.00 ▇▁▁▁▁
BoilGravity 0 1 33.31 13.50 1.00 26.00 31.50 38.00 119.00 ▂▇▁▁▁
FG 0 1 1.01 0.00 1.00 1.01 1.01 1.02 1.03 ▂▇▃▁▁
Efficiency 0 1 71.04 6.10 50.00 70.00 70.00 75.00 90.00 ▁▂▇▅▁
IBU 0 1 52.98 26.62 0.00 33.86 46.78 68.24 149.47 ▃▇▅▁▁
ABV 0 1 6.18 1.22 3.66 5.28 6.04 6.84 11.51 ▃▇▃▁▁

9 Discretización de datos

# ¿Con qué variables tendría sentido un proceso de discretización?
apply(top_top,2, function(x) length(unique(x)))
##       Clase       Color BoilGravity          FG  Efficiency         IBU 
##           5         666          87          29          41         869 
##         ABV 
##         414
# Separamos el objetivo con el resto de objetivo. El objetivo es la clase de cerveza 
top.resto <- top_top[,c(2:7)]

9.1 Normalización datos Top

normalize <- function(x){
  return ((x-min(x))/(max(x)-min(x)))
}

top.resto$Color  <- normalize(top.resto$Color)
top.resto$BoilGravity <- normalize(top.resto$BoilGravity)
top.resto$FG <- normalize(top.resto$FG)
top.resto$Efficiency <- normalize(top.resto$Efficiency)
top.resto$IBU <- normalize(top.resto$IBU)
top.resto$ABV <- normalize(top.resto$ABV)

10 Matriz de Dispersión, Histograma y Correlación de las variables seleccionadas y normalizadas.

ggpairs(top.resto, title ="Correlograma con ggpairs")

# Nice visualization of correlations
ggcorr(top.resto, method = c("everything", "pearson"))

# Scatterplot de Variables

textscatter <- function(data, mapping, ...) {
   ggplot(data, mapping, ...) + geom_text()
}

ggpairs(
  top_top, 
  title="Scatterplot de Variables",
  columns = c(2:7),
  mapping=ggplot2::aes(colour = Clase, label = ABV))

  lower = list(continuous = textscatter)

Scatterplor FG/BoilGravity

textscatter <- function(data, mapping, ...) {
   ggplot(data, mapping, ...) + geom_text()
}

ggpairs(
  top_top, 
  title="Scatterplot of Variables",
  columns = c(3:4),
  mapping=ggplot2::aes(colour = Clase, label = ABV))

  lower = list(continuous = textscatter)

Scatterplor ABV/IBU

textscatter <- function(data, mapping, ...) {
   ggplot(data, mapping, ...) + geom_text()
}

ggpairs(
  top_top, 
  title="Scatterplot of Variables",
  columns = c(6,7),
  mapping=ggplot2::aes(colour = Clase, label = ABV))

  lower = list(continuous = textscatter)
rm(correl, correl_df, correl_df_a)

Matriz Dispersión, Histograma y Correlación de datos límipios

pairs.panels(top.resto, pch=21,main="Matriz de Dispersión, Histograma y Correlación")

# Seleccionamos los de mayor corrolación
pairs.panels(top.resto[,c("BoilGravity","FG")], pch=21,main="Matriz de Dispersión, Histograma y Correlación - Boilgravity/FG")

pairs.panels(top.resto[,c("IBU","ABV")], pch=21,main="Matriz de Dispersión, Histograma y Correlación - IBU/ABV")

11 PCA

Aplicaremos PCA a las variables continuas y usaremos la variable categórica para visualizar los PCs más tarde. Noten que en el siguiente código aplicamos una transformación de registro a las variables continuas como sugiere [1] y establecemos el centro y la escala. igual a VERDADERO en la llamada a prcomp para estandarizar las variables antes de la aplicación de PCA:

stats.pca <- prcomp(top.resto) 
names(stats.pca)
## [1] "sdev"     "rotation" "center"   "scale"    "x"

11.1 Analizamos los datos

# # Aplicar PCA - escala. = VERDADERO es altamente El método de impresión devuelve la desviación estándar de cada uno de los cuatro PCs, y su rotación (o cargas), que son los coeficientes de las combinaciones lineales de las variables continuas.
# print method
print(stats.pca)
## Standard deviations (1, .., p=6):
## [1] 0.27246224 0.21774781 0.15415478 0.11938974 0.09904738 0.08399184
## 
## Rotation (n x k) = (6 x 6):
##                     PC1         PC2         PC3        PC4         PC5
## Color       -0.92255229  0.17346261 -0.03487347  0.3291019 -0.08309952
## BoilGravity -0.19055890 -0.25768270  0.13296814 -0.2834284 -0.12049066
## FG          -0.27054103 -0.27185737  0.08001683 -0.4335778  0.79118054
## Efficiency   0.04384920 -0.03061757  0.96748654  0.2372918  0.02165595
## IBU          0.08397116 -0.71338183 -0.18384660  0.6452162  0.17124509
## ABV         -0.17440822 -0.56546693  0.06979142 -0.3882774 -0.56816895
##                     PC6
## Color       -0.04876709
## BoilGravity  0.88586498
## FG          -0.18039432
## Efficiency  -0.06582728
## IBU          0.06787383
## ABV         -0.41398408

El método de trazado devuelve un trazado de las variaciones (eje y) asociadas a las PC (eje x). La figura que figura a continuación es útil para decidir cuántas PCs se deben retener para un análisis más detallado. En este simple caso con sólo 4 PCs esto no es una tarea difícil y podemos ver que las dos primeras PCs explican la mayor parte de la variabilidad de los datos.

# plot method
plot(stats.pca$x[, 1], stats.pca$x[, 2], col = top_top$Clase, main = "PCA", xlab = "PC1", ylab = "PC2")

plot(stats.pca, type = "l")

El método de resumen describe la importancia de las PC. La primera fila describe de nuevo la desviación estándar asociada a cada PC. La segunda fila muestra la proporción de la varianza en los datos explicados por cada componente, mientras que la tercera fila describe la proporción acumulativa de la varianza explicada. Podemos ver allí que los tres primeros PCs representan más del {97%} de la varianza de los datos.

# summary method
summary(stats.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6
## Standard deviation     0.2725 0.2177 0.1542 0.11939 0.09905 0.08399
## Proportion of Variance 0.4205 0.2686 0.1346 0.08074 0.05557 0.03996
## Cumulative Proportion  0.4205 0.6891 0.8237 0.90447 0.96004 1.00000

Podemos usar la función de predicción si observamos nuevos datos y queremos predecir los valores de sus PCs. Sólo para ilustrar, imaginemos que las dos últimas filas de los datos del dataset acaban de llegar y queremos ver cuáles son los valores de sus PCs:

# Predict PCs
# Preción PC de las últimas 2 lineas
predict(stats.pca, 
        newdata=tail(top.resto, 2))
##              PC1        PC2         PC3          PC4         PC5        PC6
## 10644 -0.9064684 -0.4643252  0.01369397 -0.001264468 -0.32274953 0.09897626
## 32073  0.2077856  0.1651444 -0.01923754 -0.037109444 -0.04355764 0.02706025

Grafico 2 primeros PC, que representa 84,9%

library(devtools)
install_github("vqv/ggbiplot")
## 
##   
   checking for file ‘/tmp/Rtmp9siJHw/remotes4872332e4ff8/vqv-ggbiplot-7325e88/DESCRIPTION’ ...
  
✓  checking for file ‘/tmp/Rtmp9siJHw/remotes4872332e4ff8/vqv-ggbiplot-7325e88/DESCRIPTION’
## 
  
─  preparing ‘ggbiplot’:
##    checking DESCRIPTION meta-information ...
  
✓  checking DESCRIPTION meta-information
## 
  
─  checking for LF line-endings in source and make files and shell scripts
## 
  
─  checking for empty or unneeded directories
## ─  looking to see if a ‘data/datalist’ file should be added
## ─  building ‘ggbiplot_0.55.tar.gz’
## 
  
   
## 
library(ggbiplot)
g <- ggbiplot(stats.pca, obs.scale = 1, var.scale = 1, 
              groups = top_top$Clase, ellipse = TRUE, 
              circle = TRUE)
g <- g + scale_color_discrete(name = '')
g <- g + theme(legend.direction = 'horizontal', 
               legend.position = 'top')
print(g)

### PCA on caret package

Como mencioné antes, es posible aplicar primero una transformación Box-Cox para corregir la asimetría, centrar y escalar cada variable y luego aplicar PCA en una llamada a la función de preproceso del paquete caret.

require(caret)
trans = preProcess(top_top[,2:7], 
                   method=c("BoxCox", "center", 
                            "scale", "pca"))
PC = predict(trans, top_top[,2:7])

Por defecto, la función mantiene sólo los PC necesarios para explicar al menos el 95% de la variabilidad de los datos, pero esto se puede cambiar a través del umbral de argumento.

# Retained PCs
head(PC, 3)
##              PC1        PC2         PC3         PC4       PC5       PC6
## 7608   0.9360163 -0.8551249  0.04774747  0.78069526 0.6948588 0.0215140
## 14197  1.3561292  0.9348870 -1.12542728 -0.08015398 0.5043589 0.4718226
## 2761  -0.3161090 -2.0251992  0.33484419  0.23048471 0.0604213 1.1212918
# Loadings
trans$rotation
##                     PC1         PC2         PC3         PC4        PC5
## Color       -0.29895193 -0.68156718 -0.06441744  0.61482816  0.2527764
## BoilGravity -0.51358728  0.03165505 -0.15377421 -0.45811160  0.5424035
## FG          -0.48740502 -0.29644631 -0.06401614 -0.25680399 -0.7699636
## Efficiency  -0.04901186  0.30947242 -0.90612573  0.26568134 -0.1007660
## IBU         -0.35061739  0.54256730  0.37514905  0.52362289 -0.1387185
## ABV         -0.53287661  0.23755267  0.07940508 -0.03747406  0.1402209
##                      PC6
## Color        0.005893036
## BoilGravity  0.455534598
## FG           0.108007286
## Efficiency   0.004523298
## IBU          0.385392784
## ABV         -0.795135127

11.1.1 Visualización PCA

library("factoextra")

res.pca <- prcomp(top.resto, scale = TRUE)
fviz_pca_biplot(res.pca, col.ind = top_top$Clase,
                palette = "jco", geom = "point")

La dimensión (Dim.) 1 y 2 retuvieron cerca del 61,1% (40% + 21,1%) del total de la información contenida en el conjunto de datos. Los individuos con un perfil similar se agrupan Las variables que están positivamente correlacionadas están en el mismo lado de las parcelas. Las variables que están negativamente correlacionadas están en el lado opuesto de los gráficos.

11.2 Comparamos PCA y SDV

En primer lugar, el típico PCA de las muestras sería transponer los datos de manera que las muestras sean filas de la matriz de datos.

library(rafalib)
clase <- as.character(top_top$Clase)
group <- as.fumeric(clase)
# transponemos

# x <- t(top.resto)
# Realizamos un PCA
pc <- prcomp(top.resto)
names(pc)
## [1] "sdev"     "rotation" "center"   "scale"    "x"
plot(pc$x[, 1], pc$x[, 2], col = group, main = "PCA", xlab = "PC1", ylab = "PC2")

Este PCA equivale a realizar la SVD en los datos centrados, donde el centrado ocurre en las columnas. Podemos utilizar la función de barrido para realizar operaciones arbitrarias en las filas y columnas de una matriz. El segundo argumento especifica que queremos operar en las columnas (1 se usaría para las filas), y el tercero y cuarto argumentos especifican que queremos restar la media de las columnas.

cx <- sweep(top.resto, 2, colMeans(top.resto), "-")
sv <- svd(cx)
names(sv)
## [1] "d" "u" "v"
plot(sv$u[, 1], sv$u[, 2], col = group, main = "SVD", xlab = "U1", ylab = "U2")

Así que las columnas de U de la SVD corresponden a los componentes principales x en el PCA. Además, la matriz V de la SVD equivale a la matriz de rotación devuelta por el PCA.

plot(sv$u[, 1], sv$u[, 2], col = group, main = "SVD", xlab = "U1", ylab = "U2")

plot(pc$x[, 1], pc$x[, 2], col = group, main = "PCA", xlab = "PC1", ylab = "PC2")

plot(sv$u[, 1], sv$u[, 2], col = group, main = "SVD", xlab = "U1", ylab = "U2")

Así que las columnas de U de la SVD corresponden a los componentes principales x en el PCA. Además, la matriz V de la SVD equivale a la matriz de rotación devuelta por el PCA.

sv$v[1:5, 1:5]
##             [,1]        [,2]        [,3]       [,4]        [,5]
## [1,] -0.92255229  0.17346261 -0.03487347  0.3291019 -0.08309952
## [2,] -0.19055890 -0.25768270  0.13296814 -0.2834284 -0.12049066
## [3,] -0.27054103 -0.27185737  0.08001683 -0.4335778  0.79118054
## [4,]  0.04384920 -0.03061757  0.96748654  0.2372918  0.02165595
## [5,]  0.08397116 -0.71338183 -0.18384660  0.6452162  0.17124509
pc$rotation[1:5, 1:5]
##                     PC1         PC2         PC3        PC4         PC5
## Color       -0.92255229  0.17346261 -0.03487347  0.3291019 -0.08309952
## BoilGravity -0.19055890 -0.25768270  0.13296814 -0.2834284 -0.12049066
## FG          -0.27054103 -0.27185737  0.08001683 -0.4335778  0.79118054
## Efficiency   0.04384920 -0.03061757  0.96748654  0.2372918  0.02165595
## IBU          0.08397116 -0.71338183 -0.18384660  0.6452162  0.17124509

Los elementos diagonales de D de la SVD son proporcionales a las desviaciones estándar devueltas por PCA. La diferencia es que las desviaciones estándar del prcomp son desviaciones estándar de la muestra (el prcomp devuelve estimaciones no sesgadas de la varianza de la muestra, por lo tanto con la corrección n/(n-1)). Los elementos de D se forman tomando la suma de los cuadrados de los componentes principales pero sin dividirlos por el tamaño de la muestra.

print("SDV")
## [1] "SDV"
head(sv$d^2)
## [1] 74.161437 47.366695 23.739932 14.239655  9.800573  7.047574
print("PCA")
## [1] "PCA"
head(pc$sdev^2)
## [1] 0.074235673 0.047414109 0.023763696 0.014253909 0.009810383 0.007054629
head(sv$d^2/(ncol(top.resto) - 1))
## [1] 14.832287  9.473339  4.747986  2.847931  1.960115  1.409515

Dividiendo las variaciones por la suma, obtenemos un gráfico de la relación de variación explicada por cada componente principal.

plot(sv$d^2/sum(sv$d^2), xlim = c(0, 7), type = "b", pch = 16, xlab = "Componentes principales", 
    ylab = "Varianza explicada")

plot(sv$d^2/sum(sv$d^2), type = "b", pch = 16, xlab = "Componentes principales", 
    ylab = "Varianza explicada")

Obsérvese que, al no centrar los datos antes de ejecutar el svd, se obtiene un trazado ligeramente diferente:

svNoCenter <- svd(top.resto)
plot(pc$x[, 1], pc$x[, 2], col = group, main = "PCA", xlab = "PC1", ylab = "PC2")
points(0, 0, pch = 3, cex = 4, lwd = 4)

plot(svNoCenter$u[, 1], svNoCenter$u[, 2], col = group, main = "SVD no centrada", 
    xlab = "U1", ylab = "U2")

12 Aprendizaje no supervisado

12.1 K-means

12.1.1 K-means con 4 centroides

result<- kmeans(top.resto,4) #aplico k-means con n.. de centroides(k)=4
result$size # Obtenemos el nº de grabaciones por cada cluster
## [1] 482 334  60 124
result$centers # obtenemos el valor del centro del cluster(3 centers for k=4)
##        Color BoilGravity        FG Efficiency       IBU       ABV
## 1 0.09684321   0.2247873 0.2866417  0.5266338 0.2630878 0.2380659
## 2 0.12481362   0.3116563 0.3952096  0.5324850 0.5413723 0.4333384
## 3 0.85453054   0.4862288 0.5956989  0.4962500 0.3110390 0.5362208
## 4 0.67371559   0.2594314 0.3694069  0.5211694 0.2272233 0.2354017
# Da el vector del cúmulo mostrando el cúmulo donde cae cada registro
head(result$cluster, 100)
##  7608 14197  2761 29583  1880 17415 32586  9943 23649 15226 20338  3071 21980 
##     1     1     4     2     4     1     4     1     1     4     4     1     1 
## 25647 28886 32395 27374 15655 17798  3795 28675  5658 24780 13935  2355 25867 
##     4     1     1     1     2     1     1     2     4     1     2     2     2 
## 17291  2025  4089  9282 11598 11841 25085 22616 34038 16016 13827 30786 31403 
##     1     2     2     3     1     2     4     2     2     1     4     3     1 
## 16640  7310 11614  3626  3630 21340  2475 31262 12355 23463 14871 32235 33416 
##     2     1     2     4     4     3     1     2     1     2     2     3     1 
## 34243 11665 19445 19110 31182 11083 13302 28349 29184  7381 22866 14257 30127 
##     1     1     1     4     1     4     1     1     1     2     1     4     1 
##  6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338 
##     2     2     1     2     1     1     1     4     4     2     1     1     1 
## 34258   193  9535 11829 34074 24806  9389 32397 26466   138 32918  4230 32518 
##     1     2     2     1     1     1     2     3     2     2     1     3     1 
## 30764  9581 31086 27054    53 12057 14422 12010  3190 
##     2     1     1     2     2     1     2     2     2

Verificamos el resultado del clustering

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,2)], col=result$cluster)
plot(top.resto[c(1,2)], col=top_top$Clase)
plot(top.resto[c(3,4)], col=result$cluster)
plot(top.resto[c(3,4)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(5,6)], col=result$cluster)
plot(top.resto[c(5,6)], col=top_top$Clase)
plot(top.resto[c(1,6)], col=result$cluster)
plot(top.resto[c(1,6)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,3)], col=result$cluster)
plot(top.resto[c(1,3)], col=top_top$Clase)
plot(top.resto[c(2,4)], col=result$cluster)
plot(top.resto[c(2,4)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,5)], col=result$cluster)
plot(top.resto[c(1,5)], col=top_top$Clase)
plot(top.resto[c(3,5)], col=result$cluster)
plot(top.resto[c(3,5)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,4)], col=result$cluster)
plot(top.resto[c(1,4)], col=top_top$Clase)
plot(top.resto[c(3,2)], col=result$cluster)
plot(top.resto[c(3,2)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(6,2)], col=result$cluster)
plot(top.resto[c(6,2)], col=top_top$Clase)
plot(top.resto[c(5,4)], col=result$cluster)
plot(top.resto[c(5,4)], col=top_top$Clase)

table(result$cluster,top_top$Clase)
##    
##     ALE IPA PALE PORTER STOUT
##   1  76 169  234      3     0
##   2   0 326    7      0     1
##   3   0   0    0     10    50
##   4   0   0    0     53    71

Según la tabla, el cluster 1 estan IPA, en la 2, PORTER y STOUT en la 3 ninguna, y en la cuarta ALE y PALE

El número total de instancias correctamente clasificadas son: 786 El número total de casos clasificados incorrectamente son: 214 Exactitud = 786/(1000) = 0.786 i.e. nuestro modelo ha alcanzado un 78,6% de exactitud!

12.1.2 K-means con 3 centroides

result<- kmeans(top.resto,3) #aplico k-means con n.. de centroides(k)=3
result$size # Obtenemos el nº de grabaciones por cada cluster
## [1]  60 126 814
result$centers # obtenemos el valor del centro del cluster(3 centers for k=3)
##       Color BoilGravity        FG Efficiency       IBU       ABV
## 1 0.8545305   0.4862288 0.5956989  0.4962500 0.3110390 0.5362208
## 2 0.6687941   0.2569276 0.3684076  0.5228175 0.2255414 0.2336467
## 3 0.1076644   0.2607338 0.3311405  0.5287930 0.3776218 0.3184682
# Da el vector del cúmulo mostrando el cúmulo donde cae cada registro
head(result$cluster, 100)
##  7608 14197  2761 29583  1880 17415 32586  9943 23649 15226 20338  3071 21980 
##     3     3     2     3     2     3     2     3     3     2     2     3     3 
## 25647 28886 32395 27374 15655 17798  3795 28675  5658 24780 13935  2355 25867 
##     2     3     3     3     3     3     3     3     2     3     3     3     3 
## 17291  2025  4089  9282 11598 11841 25085 22616 34038 16016 13827 30786 31403 
##     3     3     3     1     3     3     2     3     3     3     2     1     3 
## 16640  7310 11614  3626  3630 21340  2475 31262 12355 23463 14871 32235 33416 
##     3     3     3     2     2     1     3     3     3     3     3     1     3 
## 34243 11665 19445 19110 31182 11083 13302 28349 29184  7381 22866 14257 30127 
##     3     3     3     2     3     2     3     3     3     3     3     2     3 
##  6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338 
##     3     3     3     3     3     3     3     2     2     3     3     3     3 
## 34258   193  9535 11829 34074 24806  9389 32397 26466   138 32918  4230 32518 
##     3     3     3     3     3     3     3     1     3     3     3     1     3 
## 30764  9581 31086 27054    53 12057 14422 12010  3190 
##     3     3     3     3     3     3     3     3     3

Verificamos el resultado del clustering

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,2)], col=result$cluster)
plot(top.resto[c(1,2)], col=top_top$Clase)
plot(top.resto[c(3,4)], col=result$cluster)
plot(top.resto[c(3,4)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(5,6)], col=result$cluster)
plot(top.resto[c(5,6)], col=top_top$Clase)
plot(top.resto[c(1,6)], col=result$cluster)
plot(top.resto[c(1,6)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,3)], col=result$cluster)
plot(top.resto[c(1,3)], col=top_top$Clase)
plot(top.resto[c(2,4)], col=result$cluster)
plot(top.resto[c(2,4)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,5)], col=result$cluster)
plot(top.resto[c(1,5)], col=top_top$Clase)
plot(top.resto[c(3,5)], col=result$cluster)
plot(top.resto[c(3,5)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,4)], col=result$cluster)
plot(top.resto[c(1,4)], col=top_top$Clase)
plot(top.resto[c(3,2)], col=result$cluster)
plot(top.resto[c(3,2)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(6,2)], col=result$cluster)
plot(top.resto[c(6,2)], col=top_top$Clase)
plot(top.resto[c(5,4)], col=result$cluster)
plot(top.resto[c(5,4)], col=top_top$Clase)

table(result$cluster,top_top$Clase)
##    
##     ALE IPA PALE PORTER STOUT
##   1   0   0    0     10    50
##   2   0   0    0     55    71
##   3  76 495  241      1     1

Según la tabla, el cluster 1 estan ALE y PALE, en la 2, claramente IPA y en la tercera, PORTER y STOUT.

El número total de instancias correctamente clasificadas son: 835 El número total de casos clasificados incorrectamente son: 165 Exactitud = 835/(835+165) = 0.835 i.e. nuestro modelo ha alcanzado un 83,5% de exactitud!

12.1.3 K-means con 2 centroides

result<- kmeans(top.resto,2) #aplico k-means con n.. de centroides(k)=2
result$size # Obtenemos el nº de grabaciones por cada cluster
## [1] 183 817
result$centers # obtenemos el valor del centro del cluster(3 centers for k=2)
##       Color BoilGravity        FG Efficiency       IBU       ABV
## 1 0.7343862   0.3340048 0.4438569  0.5131148 0.2545376 0.3342940
## 2 0.1086733   0.2602950 0.3310696  0.5289933 0.3768474 0.3178336
# Da el vector del cúmulo mostrando el cúmulo donde cae cada registro
head(result$cluster, 100)
##  7608 14197  2761 29583  1880 17415 32586  9943 23649 15226 20338  3071 21980 
##     2     2     1     2     1     2     1     2     2     1     1     2     2 
## 25647 28886 32395 27374 15655 17798  3795 28675  5658 24780 13935  2355 25867 
##     1     2     2     2     2     2     2     2     1     2     2     2     2 
## 17291  2025  4089  9282 11598 11841 25085 22616 34038 16016 13827 30786 31403 
##     2     2     2     1     2     2     1     2     2     2     1     1     2 
## 16640  7310 11614  3626  3630 21340  2475 31262 12355 23463 14871 32235 33416 
##     2     2     2     1     1     1     2     2     2     2     2     1     2 
## 34243 11665 19445 19110 31182 11083 13302 28349 29184  7381 22866 14257 30127 
##     2     2     2     1     2     1     2     2     2     2     2     1     2 
##  6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338 
##     2     2     2     2     2     2     2     1     1     2     2     2     2 
## 34258   193  9535 11829 34074 24806  9389 32397 26466   138 32918  4230 32518 
##     2     2     2     2     2     2     2     1     2     2     2     1     2 
## 30764  9581 31086 27054    53 12057 14422 12010  3190 
##     2     2     2     2     2     2     2     2     2

Verificamos el resultado del clustering

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,2)], col=result$cluster)
plot(top.resto[c(1,2)], col=top_top$Clase)
plot(top.resto[c(3,4)], col=result$cluster)
plot(top.resto[c(3,4)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(5,6)], col=result$cluster)
plot(top.resto[c(5,6)], col=top_top$Clase)
plot(top.resto[c(1,6)], col=result$cluster)
plot(top.resto[c(1,6)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,3)], col=result$cluster)
plot(top.resto[c(1,3)], col=top_top$Clase)
plot(top.resto[c(2,4)], col=result$cluster)
plot(top.resto[c(2,4)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,5)], col=result$cluster)
plot(top.resto[c(1,5)], col=top_top$Clase)
plot(top.resto[c(3,5)], col=result$cluster)
plot(top.resto[c(3,5)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(1,4)], col=result$cluster)
plot(top.resto[c(1,4)], col=top_top$Clase)
plot(top.resto[c(3,2)], col=result$cluster)
plot(top.resto[c(3,2)], col=top_top$Clase)

par(mfrow=c(2,2), mar=c(5,4,2,2))
plot(top.resto[c(6,2)], col=result$cluster)
plot(top.resto[c(6,2)], col=top_top$Clase)
plot(top.resto[c(5,4)], col=result$cluster)
plot(top.resto[c(5,4)], col=top_top$Clase)

table(result$cluster,top_top$Clase)
##    
##     ALE IPA PALE PORTER STOUT
##   1   0   0    0     63   120
##   2  76 495  241      3     2

Según la tabla, el cluster 2 estan ALE, IPA Y PALE y el cluster 1 están PORTER Y STOUT

El número total de instancias correctamente clasificadas son: 995 El número total de casos clasificados incorrectamente son: 5 Exactitud = 995/(1000) = 0.995 i.e. nuestro modelo ha alcanzado un 99,% de exactitud!

12.2 Busqueda del nº de clusters

library(factoextra)
fviz_nbclust(x = top.resto, FUNcluster = kmeans, method = "wss", k.max = 7, 
             diss = get_dist(top.resto, method = "euclidean"), nstart = 30)

# Este mismo análisis también puede realizarse sin recurrir a la función fviz_nbclust().

calcular_totwithinss <- function(n_clusters, datos, iter.max=1000, nstart=50){
  # Esta función aplica el algoritmo kmeans y devuelve la suma total de
  # cuadrados internos.
  cluster_kmeans <- kmeans(centers = n_clusters, x = top.resto, iter.max = iter.max,
                           nstart = nstart)
  return(cluster_kmeans$tot.withinss)
}
# library(purrr)
# Se aplica esta función con para diferentes valores de k
total_withinss <- map_dbl(.x = 1:7,
                          .f = calcular_totwithinss,
                          datos = top.resto)
total_withinss
## [1] 176.35587 112.79074  85.87822  76.58989  69.04624  64.00273  60.15495
data.frame(n_clusters = 1:7, suma_cuadrados_internos = total_withinss) %>%
  ggplot(aes(x = n_clusters, y = suma_cuadrados_internos)) +
    geom_line() +
    geom_point() +
    scale_x_continuous(breaks = 1:7) +
    labs(title = "Evolución de la suma total de cuadrados intra-cluster") +
    theme_bw()

12.3 PAM

library(cluster)
library(factoextra)
fviz_nbclust(x = top.resto, FUNcluster = pam, method = "wss", k.max = 5,
             diss = dist(top.resto, method = "manhattan"))

set.seed(123)
pam_clusters <- pam(x = top.resto, k = 3, metric = "manhattan")
pam_clusters
## Medoids:
##        ID      Color BoilGravity        FG Efficiency       IBU       ABV
## 26633 384 0.04732984   0.2076271 0.2903226        0.5 0.2466716 0.2127389
## 7922  456 0.79057592   0.3135593 0.4193548        0.5 0.2511541 0.3617834
## 8841  727 0.10198953   0.3050847 0.3870968        0.5 0.4388841 0.4433121
## Clustering vector:
##  7608 14197  2761 29583  1880 17415 32586  9943 23649 15226 20338  3071 21980 
##     1     1     2     3     2     1     2     3     3     2     2     1     3 
## 25647 28886 32395 27374 15655 17798  3795 28675  5658 24780 13935  2355 25867 
##     2     1     1     1     3     1     1     3     1     1     3     3     3 
## 17291  2025  4089  9282 11598 11841 25085 22616 34038 16016 13827 30786 31403 
##     1     3     3     2     1     3     2     3     3     1     2     2     1 
## 16640  7310 11614  3626  3630 21340  2475 31262 12355 23463 14871 32235 33416 
##     3     1     3     2     2     2     1     3     1     3     3     2     1 
## 34243 11665 19445 19110 31182 11083 13302 28349 29184  7381 22866 14257 30127 
##     1     1     1     1     1     2     1     1     3     3     3     2     1 
##  6822 21654 29721 25478 10901 31898 21894 23807 21576 22799 26038 23601 14338 
##     3     3     1     3     1     1     1     1     2     3     1     3     3 
## 34258   193  9535 11829 34074 24806  9389 32397 26466   138 32918  4230 32518 
##     1     3     3     1     1     1     1     2     3     3     1     2     1 
## 30764  9581 31086 27054    53 12057 14422 12010  3190 32889 28666 33777  5927 
##     3     3     1     3     3     1     3     3     3     2     3     1     1 
##  6298 21656  3586 28856 31152 19089 13566 33631 29572 12154 26093 25175 26830 
##     3     1     3     3     3     3     3     1     3     3     3     3     1 
## 11433  4735 23258   734  4014 28804 18054 12340 34510 12765 15247 20063  4216 
##     1     1     3     3     2     3     1     2     2     1     1     3     3 
## 12138 14321 31450 31817 21195  1331 25300 29627 17785 24632  9379 17188 18878 
##     1     3     1     3     2     3     1     3     1     1     1     3     2 
## 11621 20658 11331  5421 27666 31845   672 34152  6760 24223 26155  2771 24749 
##     3     1     1     3     2     2     3     1     1     1     3     1     1 
## 23819 32661  6924  9848   241 19258  2578 28316  5993 29561 10786 27129  9976 
##     1     3     1     3     3     3     3     1     3     2     3     3     2 
## 18389 12964  2743 20800 12835  4775  4497 31138 16385 10368 17751 13497 33188 
##     3     1     2     1     2     1     2     2     3     1     1     3     1 
##   449 31363 33331 20086 34142 25046 23742  3112  7969 14044 27219 21428  8305 
##     3     3     3     3     3     1     1     3     3     3     1     2     3 
## 30962 18861  3222 27706 16803 12674 10451 30355 12300   749 26959 22037 22332 
##     1     1     1     1     2     2     2     3     3     1     2     2     3 
## 24527 11421 32385 21988  8586 24383  3988  6144 32093 25432 30705 23910 24156 
##     3     2     1     2     3     3     3     3     3     3     2     3     1 
## 30235 12158  5119 22779 21311  5409 31292 23127 14005  7136 32610 33622  5753 
##     1     2     1     1     1     3     2     3     1     1     3     3     1 
## 11095  2219 24759  2167 33920  5987 27649 13207 15431  9078  4292 29265 18837 
##     3     2     3     3     1     3     3     2     1     1     1     2     1 
## 34205 16154  5650  2359 20160 11844  4684  2413  7702  7753 25268 33966 20429 
##     2     1     3     2     1     3     1     1     1     3     1     3     1 
##  6966  3784 11101 31514 29685 20459  3996 13748  2712   799  5682 17454 33765 
##     3     1     1     1     1     1     1     1     3     2     3     3     3 
##  2927 30444  6849 24469 10827 23945 23094 30561 31923 16407 32378  8860 27048 
##     1     1     1     2     1     1     3     1     1     1     1     1     3 
## 12301 31516 32134 14358 21520 16575 28474 18829  6881 18428 29492 23394 22384 
##     2     2     1     3     1     1     3     3     1     1     2     1     2 
## 16566 34144  9663 22290 23107 17819  9574 31926  4287 17859  3117 15663 21102 
##     1     2     2     3     3     3     2     3     3     2     1     1     3 
##  8852 34361  2239  9881 32863 13253 23360  9605 28504  6177 28472  2739 26559 
##     2     3     1     3     1     3     3     1     3     3     3     3     3 
## 34249 29130 14908  7097 32348  7342  7168  5087 24425   912  4856 10280 16888 
##     1     3     3     1     1     3     3     2     3     3     1     3     1 
##  3673 28685 33271  2857  4669 10363 17876 11955   344 27611 10323 33720 12092 
##     1     2     1     2     2     1     3     3     3     1     3     1     2 
##  6961 31781 25176 13739 27654 33210 16625  6860 20316 12988 18670   229 32474 
##     1     3     1     3     2     3     3     3     3     2     1     1     1 
## 15745 17705 31522 27738  9075  8696  8437 32364 12408 30125 25362 25771 10373 
##     3     1     1     3     2     3     2     1     3     3     1     3     3 
## 18151 19147 25206  1586 33320  6839 26633 27185  1982 31217  2546 25673  6753 
##     1     1     1     3     2     1     1     1     1     1     1     1     2 
##  1516 33899 21318  1650 22326 16735 34328 18300 31964 18269 34259 34500 23192 
##     1     1     2     3     1     3     1     2     1     1     1     1     1 
## 18554 18599  8941 31283 20189 10184 16456 16760 19687 28208 17521  9814 14886 
##     1     1     3     2     1     1     1     1     3     3     2     1     2 
## 12202  3524 18456 10576 30222 32458  8788  3771 13870 28041 15369 21909 26111 
##     2     3     3     3     3     3     1     1     3     2     1     3     3 
## 13012 24199 20178 18819 20663 24672 18284 10767 16083  6724 16059  2213 31050 
##     2     1     1     1     1     3     3     3     1     2     3     3     2 
## 20518 27730  3258 12779  3893 17437  8611   808 32463  5017 33636  4679 33911 
##     3     3     2     1     1     2     2     2     3     1     3     1     3 
##  7922 10060  6767 16630 32680 26584 27149 17364 26270  5337 25240 11554 23935 
##     2     1     3     1     1     3     1     3     1     1     1     3     3 
## 26947  9753 25245 23613 20092 30595 34651 20098 25328  8088 13914 26013   113 
##     1     3     1     3     1     1     2     1     3     3     1     3     1 
##  9368 25151 31293 25800 10883 30560 28006 10770 20439 27880 28945 26846  7910 
##     1     3     1     1     1     1     3     3     1     3     3     1     1 
##  7971  9707  2941 27802 18831 14247 26932 18590 27759  1106  8285 14588  5787 
##     3     1     2     3     3     1     3     2     1     1     2     3     2 
## 11822 19806  8518 31755 21095  6117  9894 25544  2577  2937 23963 25230 21984 
##     2     2     2     3     3     3     2     1     1     3     1     1     1 
## 22917 15149   648 26141 12425 11320 24764 25540 22791 27367  1934 10054 20828 
##     1     1     1     1     1     1     3     1     1     1     1     1     1 
## 24978 20561  4180 30896 34345 32503  1512 11736 16700  9222 33104 23793  5656 
##     3     3     3     3     3     2     1     3     3     2     3     3     1 
## 19978 12857  2586 31307 22622  9926 33146 19340  9327 20583 34076 17956  3155 
##     3     1     3     3     1     3     3     1     3     1     1     3     1 
## 11329 19165 27103  8261 27224 26036 24483 19154 30287 22044 31647 11449  5854 
##     3     3     2     3     1     3     2     1     1     1     2     1     3 
## 29166 13391 26777 16044 13230  3601 12120 19612  6920 13470 16419 29386 21247 
##     3     2     1     3     1     3     3     1     3     1     1     3     1 
## 22963  2407 20824 13588 30858 20977 10299 31482  7709 11260 15679 23758 11819 
##     1     1     1     1     3     3     2     2     1     2     1     1     3 
## 25045 25691 34045 14862  6018 25379  1631  9545 17387  6173 23761 22724  6163 
##     1     3     1     1     1     1     3     3     3     1     2     1     1 
## 12446 20853 24424 23673 19491 23422 11028  8215 28104  8861 18763 30570 17361 
##     1     3     1     1     3     3     3     1     1     3     3     3     2 
## 32374 23374 20132 32399 15878   445 16729  5732   412  1499 17235 14100 25947 
##     1     1     3     3     1     2     1     2     2     3     1     1     2 
## 10292  7000  5629 16379 22638 13612 20462 25533   570 29024  3919 30728 32787 
##     1     3     1     3     3     3     1     1     3     3     1     2     1 
## 23410 27516 25380 12502 25604 18898  3324  8595 19931  7447 23501 24258  2042 
##     3     1     3     1     1     1     2     3     1     2     1     3     3 
##  3682 22571 11644 26723 22669  9939 27331 21300 23171 15569 27262 28941 16608 
##     1     3     3     3     3     3     1     1     1     3     3     1     3 
##  9562 16025 32026 21892 19136 30941  3664  3678 11199 11733  2095 32237 28360 
##     3     2     3     1     3     1     2     3     1     3     2     1     2 
## 32638 17266  2119 14224 21716  6504  7012 34431 17181  5666  5932 31973 20534 
##     1     2     3     1     3     3     3     1     2     2     1     1     2 
## 23046 32488 21125   931 33409 11610 25964 20894 15437  1391  6637 32265 27052 
##     3     3     3     3     1     1     1     1     1     1     2     3     3 
## 26027 10040 26016  4510 34448 21468  7739 21162  3162 19627 26508  8841 21940 
##     1     1     1     3     3     2     3     1     3     3     3     3     3 
## 16547   513 25489 17103 11429 19030 11810 21505  9998 29860  9452 17877 27508 
##     3     1     1     3     1     1     1     3     1     1     1     3     1 
## 21828 12626  5619  3097 16031  6711  4447  1036 34098 24406 10664 21200 34445 
##     1     1     2     3     1     3     1     3     1     3     1     1     1 
## 32075 32165 23704 24766  3070 34217  8925 28847 28553 15118 34194  4695  5513 
##     2     1     3     1     1     1     3     2     3     3     2     1     1 
## 13166 12588  9357 34491 23676 25283 15864 11864  7275  2465 31931 29533 19697 
##     3     3     3     1     2     3     2     3     3     1     1     1     1 
## 25759 16995 17366 30472 22184 12955 16064  1695 29364  6028   916  7557 32333 
##     2     1     1     3     3     3     1     3     3     1     3     3     3 
## 31783 31452 21752 20120 23125 33530  1110 25730 15059 19844  1553 32371 17588 
##     3     2     2     1     2     1     2     1     3     1     3     1     1 
##  4538 32202 14249 17447 25373  5766   911   526 27038 14905 28382 28460 27889 
##     3     1     2     2     1     3     3     3     1     1     3     3     2 
## 11232  4830  2874 11386  2582 34027 14743 29493 20137  7468 30676 12264 23159 
##     1     1     3     3     3     2     1     2     1     1     3     2     1 
## 15354 24828 29744 27474 24445 27679  9669 24582 10418 10283  7766  3634 12453 
##     3     3     1     1     1     2     1     1     3     1     3     1     2 
## 30511 33908 18801 14672 30276 15849 13151 16768 17487  5001 26145     3 19837 
##     1     1     1     3     1     3     2     2     1     1     1     1     3 
## 15914 14991  4345 30800 27642 24956 20970 12007  7906 11214 30295 15281 17436 
##     1     3     3     1     3     2     1     3     3     1     1     1     3 
## 31227 27582 10670 32496  5663 17241  3780 11357 34009 19011  5978  8703   823 
##     3     3     2     1     3     2     1     1     3     3     3     1     2 
## 14090 10016  3186 12099 17328 10395 15366 17275 28546 27587   125  5287 11253 
##     3     3     3     1     2     3     1     1     3     1     1     1     1 
## 26709   281 31011 32774 32554 34455 25886  3160  1733 24047  6794 15927  7337 
##     3     1     3     1     2     3     2     3     3     1     3     2     1 
## 14665 33114 14463  8330 30046 31651 34399  6738 17823  9090 15352  2698   552 
##     1     3     1     1     1     3     1     1     3     3     1     1     3 
## 18355 17726 32563  1155 20706 33276 26471 11263 14454 22573 27896 29512 28705 
##     3     3     2     2     3     1     2     1     1     1     1     1     3 
## 32753 10429  6939 34636  1671 16005 11891  2984 12058 19912 13816  8441 19128 
##     1     2     3     3     3     2     1     3     1     3     3     2     3 
## 13181 30698 31294 10625 14264 22777 20290 23830 30499 31949  5488 22444 17920 
##     3     1     1     3     1     3     1     1     3     3     1     3     3 
## 14478  5895 15246  8264  8056 22496 14190 26273 31471  9079  5768 31649 18404 
##     3     2     3     3     3     3     3     1     1     1     1     2     3 
## 20520 32982 16987 18360 33832  7999  1628 30171 11393  5456  7231 13510 26937 
##     3     1     1     2     3     3     1     1     1     3     1     1     1 
##  9155 21569 13945 10662 17172 14920 13940 11413 34070 10387 10644 32073 
##     2     3     1     3     1     1     1     3     3     2     2     1 
## Objective function:
##     build      swap 
## 0.5337826 0.5104537 
## 
## Available components:
##  [1] "medoids"    "id.med"     "clustering" "objective"  "isolation" 
##  [6] "clusinfo"   "silinfo"    "diss"       "call"       "data"
fviz_cluster(object = pam_clusters, data = top.resto, ellipse.type = "t",
             repel = TRUE) +
  theme_bw() +
  labs(title = "Resultados clustering PAM") +
  theme(legend.position = "none")

## CLARA

library(cluster)
library(factoextra)
clara_clusters <- clara(x = top.resto, k = 3, metric = "manhattan", stand = TRUE,
                        samples = 10, pamLike = TRUE)
clara_clusters
## Call:     clara(x = top.resto, k = 3, metric = "manhattan", stand = TRUE,      samples = 10, pamLike = TRUE) 
## Medoids:
##            Color BoilGravity        FG Efficiency       IBU       ABV
## 32073 0.05089005   0.2288136 0.2258065      0.500 0.2280725 0.2178344
## 12264 0.65005236   0.2500000 0.4032258      0.500 0.2249281 0.3006369
## 31227 0.07790576   0.2966102 0.3548387      0.575 0.4129926 0.3324841
## Objective function:   4.546298
## Clustering vector:    Named int [1:1000] 1 1 2 3 2 1 2 3 3 2 2 1 3 2 1 1 3 3 ...
##  - attr(*, "names")= chr [1:1000] "7608" "14197" "2761" "29583" "1880" "17415" "32586" ...
## Cluster sizes:            352 176 472 
## Best sample:
##  [1] 3071  13935 25085 25478 9535  12138 20658 672   12964 17751 12674 6177 
## [13] 34249 33271 12408 25673 18269 28208 13012 33911 23613 20092 25151 26846
## [25] 9894  33146 13588 17387 32399 32787 3324  3682  21828 34445 34217 34194
## [37] 12955 31783 12264 30511 31227 5663  12058 13816 11413 32073
## 
## Available components:
##  [1] "sample"     "medoids"    "i.med"      "clustering" "objective" 
##  [6] "clusinfo"   "diss"       "call"       "silinfo"    "data"
fviz_cluster(object = clara_clusters, ellipse.type = "t", geom = "point",
             pointsize = 2.5) +
  theme_bw() +
  labs(title = "Resultados clustering CLARA") +
  theme(legend.position = "none")

## Hierarchical K-means clusteriong

library(factoextra)
# Se obtiene el dendrograma de hierarchical clustering para elegir el número de
# clusters.
set.seed(101)
hc_euclidea_completo <- hclust(d = dist(x = top.resto, method = "euclidean"),
                               method = "complete")
fviz_dend(x = hc_euclidea_completo, cex = 0.5, main = "Linkage completo",
          sub = "Distancia euclídea") +
  theme(plot.title =  element_text(hjust = 0.5, size = 15))

13 Referencias:

https://www.r-bloggers.com/experimentation-with-unsupervised-learning/

https://www.r-bloggers.com/computing-and-visualizing-pca-in-r/

https://askubuntu.com/questions/161778/how-do-i-use-wget-curl-to-download-from-a-site-i-am-logged-into

https://genomicsclass.github.io/book/pages/pca_svd.html

https://rpubs.com/Nitika/kmeans_Iris

https://rpubs.com/Joaquin_AR/310338

https://www.adictosaltrabajo.com/2019/11/28/deteccion-y-reemplazo-de-outliers-con-r/